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Abstract 

The reversible Markov chains that drive the data augmentation (DA) and sandwich algo- 
rithms define self-adjoint operators whose spectra encode the convergence properties of the 
algorithms. When the target distribution has uncountable support, as is nearly always the case 
in practice, it is generally quite difficult to get a handle on these spectra. We show that, if the 
augmentation space is finite, then (under regularity conditions) the operators defined by the DA 
and sandwich chains are compact, and the spectra are finite subsets of [0, 1). Moreover, we 
prove that the spectrum of the sandwich operator dominates the spectrum of the DA operator in 
the sense that the ordered elements of the former are all less than or equal to the corresponding 
elements of the latter. As a concrete example, we study a widely used DA algorithm for the ex- 
ploration of posterior densities associated with Bayesian mixture models ( |Diebolt and Robert] 



1994). In particular, we compare this mixture DA algorithm with an alternative algorithm pro- 



posed by Friihwirth-Schnatter (2001 1 that is based on random label switching. 



1 Introduction 

Suppose that fx '■ K p — > [0, oo) is a probability density function that is intractable in the sense that 
expectations with respect to fx cannot be computed analytically. If direct simulation from fx is 
infeasible, then classical Monte Carlo methods cannot be used to explore fx and one might resort 
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to a Markov chain Monte Carlo (MCMC) method such as the data augmentation (DA) algorithm 
( |Hobert| |2011[ |Liu, Wong and Kong[ |1994[ |Tanner and Wong[ |1987| ). To build a DA algorithm, 
one must identify a joint density, say / : l p x I' 4 [0, oo), that satisfies two conditions: (i) the 
x-marginal of f(x, y) is fx, and (ii) sampling from the associated conditional densities, fx\y i'lu) 
and fy\x('\ x ), i s straightforward. (The y-coordinate may be discrete or continuous.) The first of 
the two conditions allows us to construct a Markov chain having fx as an invariant density, and 
the second ensures that we are able to simulate this chain. Indeed, let {X n }^ =0 be a Markov chain 
whose dynamics are defined (implicitly) through the following two-step procedure for moving from 
the current state, X n = x, to X n+ \. 



Iteration n + 1 of the DA Algorithm: 



1. Draw Y ~ fy\x('\ x )> and call the observed value y 

2. Draw X n+1 ~ f x \y(-\y) 



It is well known and easy to establish that the DA Markov chain is reversible with respect to fx, and 
this of course implies that fx is an invariant density ( [Liu et aL 1994| ). Consequently, if the chain 
satisfies the usual regularity conditions (see Section [2]), then we can use averages to consistently es- 
timate intractable expectations with respect to fx (Tierney 1994). The resulting MCMC algorithm 
is known as a DA algorithm for fx- (Throughout this section, fx is assumed to be a probability 
density function, but starting in Section|2]a more general version of the problem is considered.) 

When designing a DA algorithm, one is free to choose any joint density that satisfies conditions 
(i) and (ii). Obviously, different joint densities will yield different DA chains, and the goal is to find 
a joint density whose DA chain has good convergence properties. (This is formalized in Section [3] 
using x 2 -distance to stationarity.) Unfortunately, the "ideal" joint density, which yields the DA 
chain with the fastest possible rate of convergence, does not satisfy the simulation requirement. 
Indeed, consider f±(x,y) = fx(x) gy(y), where gy(y) is any density function on M. q . Since 
f±(x, y) factors, f x \y (x\y) = fx{%) and it follows that the DA chain is just an iid sequence from 
fx- Of course, this ideal DA algorithm is useless from a practical standpoint because, in order to 
simulate the chain, we must draw from fx, which is impossible. We return to this example later in 
this section. 
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It is important to keep in mind that there is no inherent interest in the joint density f(x,y). It 
is merely a tool that facilitates exploration of the target density, fx{x). This is the reason why the 
DA chain does not possess a y-coordinate. In contrast, the two-variable Gibbs sampler based on 
fx\Y{'\y) an d fy\x('\ x )' which is used to explore f(x, y), has both x and y-coordinates. So, while 
the two-step procedure described above can be used to simulate both the DA and Gibbs chains, there 
is one key difference. When simulating the DA chain, we do not keep track of the y-coordinate. 

Every reversible Markov chain defines a self-adjoint operator whose spectrum encodes the con- 
vergence properties of the chain (Diaconis, Khare and Saloff-Coste 2008; Mira and Geyer| 1999 



Rosenthal 2003 1. Let X ~ fx and consider the space of functions g such that the random variable 



g(X) has finite variance and mean zero. To be more precise, define 

L 2 (f x ) = {g:M. p ^R: [ g\x) f x (x) dx < oo and / g{x) f x (x) dx = 



Let k(x'\x) be the Markov transition density (Mtd) of the DA chain. (See Section [5] for a formal 
definition). This Mtd defines an operator, K : L\ — > h\ , that maps g(x) to 



(Kg)(x) := / g{x')k(x'\x)dx . 
Jw 



Of course, {Kg){x) is just the expected value of g(X\) given that Xq = x. Let I : L\ — > L\ denote 
the identity operator, which leaves functions unaltered, and consider the operator K — XI, where 
A G M. By definition, K — XI is invertible if, for each h 6 Lq, there exists a unique g G Lq sucn 
that {{K — XI)g){x) = (Kg)(x) — Xg(x) = h(x). The spectrum of K, which we denote by Sp(K), 
is simply the set of A such that K — XI is not invertible. Because K is defined through a DA chain, 
Sp(if) C [0, 1] (see Section[3]>. The number of elements in Sp(K) may be finite, countably infinite 
or uncountable. 

In order to understand what "good" spectra look like, consider the ideal DA algorithm in- 
troduced earlier. Let k± and K± denote the Mtd and the corresponding operator, respectively. 
In the ideal case, X n+ \ is independent of X n and has density fx- Therefore, the Mtd is just 

k±(x'\x) = fx{x') and 

(K±g)(x) = / g{x) fx(x') dx' = , 



which implies that 

((K ± -XI)g)(x) = -Xg(x) 
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It follows that K±— A/ is invertible as long as A ^ 0. Hence, the "ideal spectrum" is Sp(K±) = {0}. 
Loosely speaking, the closer Sp(K) is to {0}, the faster the DA algorithm converges ( |Diaconis et al.| 
[2008] ). 

Unfortunately, in general, there is no simple method for calculating Sp(K). Even getting a 
handle on Sp(K) is currently difficult. However, there is one situation where Sp(K) has a very 
simple structure. Let Y = {y £ R q : fy(y) > 0}, where fy(y) = L p f(x, y) dx. We show that 
when Y is a finite set, Sp(K) consists of a finite number of elements that are directly related to the 
Markov transition matrix (Mtm) of the so-called conjugate chain, which is the reversible Markov 
chain that lives on Y and makes the transition y — )■ y' with probability J RP fy\x{y'\x) fx\Y( x \y) dx. 
In particular, we prove that when |Y| = d < oo, Sp(K) consists of the point {0} together with the 
d — 1 smallest eigenvalues of the Mtm of the conjugate chain. We use this result to prove that the 
spectrum associated with a particular alternative to the DA chain is closer than Sp(K) to the ideal 
spectrum, {0}. 

DA algorithms often suffer from slow convergence, which is not surprising given the close 



connection between DA and the notoriously slow to converge EM algorithm (see e.g. van Dyk and 



Meng[ 2001 ). Over the last decade, a great deal of effort has gone into modifying the DA algorithm 



to speed convergence. See, for example, |Meng and van Dyk| ( |1999| ), |Liu and Wu| ( |1999| ), |Liu and 



Sabattij ( |2000] >, |van Dyk and Meng| pOOT) , jiPapaspil iopoulos, Roberts and Skold] ( |2007] >, |Hobert 



and Marchev ( 2008[ ) and Yu and Meng (2011 ). In this paper, we focus on the so-called sandwich 



algorithm, which is a simple alternative to the DA algorithm that often converges much faster. Let 
r(y'\y) be an auxiliary Mtd (or Mtm) that is reversible with respect to fy, and consider a new 
Markov chain, { Ar n }^Lo> tnat moves from X n = x to X n+ \ via the following three-step procedure. 



Iteration n + 1 of the Sandwich Algorithm: 

1. Draw Y ~ fy\x{-\x), and call the observed value y 

2. Draw Y' ~ r(-\y), and call the observed value y' 

3. Draw X n+l ~ f x \Y{-\y') 



A routine calculation shows that the sandwich chain remains reversible with respect to fx, so it is 
a viable alternative to the DA chain. The name "sandwich algorithm" was coined by | Yu and Meng 
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(2011 1 and is based on the fact that the extra draw from r(-\y) is sandwiched between the two steps 
of the DA algorithm. Clearly, on a per iteration basis, it is more expensive to simulate the sandwich 
chain. However, it is often possible to find an r that leads to a substantial improvement in mixing 
despite the fact that it only provides a low-dimensional (and hence inexpensive) perturbation on the 
Y space. In fact, the computational cost of drawing from r is often negligible relative to the cost 
of drawing from fy\x('\ x ) an d fx\y{'\y)- Concrete examples can be found in Meng and van Dyk 



< fT999] ), |Liu and Wu| ([1999), |van Dyk and Meng| pOOT] ), |Roy and Hob"ert| ( p007| ) and Section [5] of 
this paper. 

Let k(x'\x) denote the Mtd of the sandwich chain. Also, let K and Sp(K) denote the corre- 
sponding operator and its spectrum. The main theoretical result in this paper provides conditions 
under which Sp(K) is closer than Sp(K) to the ideal spectrum. Recall that when |Y| = d < oo, 
Sp(K) consists of the point {0} and the d — 1 smallest eigenvalues of the Mtm of the conjugate 
chain. If, in addition, r is idempotent (see Section [4] for the definition), then Sp{K) consists of 
the point {0} and the d — 1 smallest eigenvalues of a different d x d Mtm, and < A,; < Aj 
for all i € {1, 2, . . . , d — 1}, where Aj and Aj are the ith largest elements of Sp(-fiT) and Sp(K), 
respectively. So Sp(i^) dominates Sp(K) in the sense that the ordered elements of Sp(K) are uni- 
formly less than or equal to the corresponding elements of Sp(K). We conclude that the sandwich 
algorithm is closer than the DA algorithm to the gold standard of classical Monte Carlo. 

One might hope for a stronger result that quantifies the extent to which the sandwich chain is 
better than the DA chain, but such a result is impossible without further assumptions. Indeed, if we 
take the auxiliary Markov chain on Y to be the degenerate chain that is absorbed at its starting point, 
then the sandwich chain is the same as the DA chain. 

To illustrate the huge gains that are possible through the sandwich algorithm, we introduce a 
new example involving a Bayesian mixture model. Let Z\, . . . , Z m be a random sample from a 
/c-component mixture density taking the form 

k 

^pjhg^z) , (1) 
3=1 

where #i, . . . , 0& G C M. 1 , {ho(-) : 9 £ 6} is a parametric family of densities, and the pjs are 
nonnegative weights that sum to one. Of course, a Bayesian analysis requires priors for the unknown 
parameters, which are 6 = (#1, #2, • • • , 9k) T and p = (pi,P2, ■ ■ ■ ,Pk) T - In typical applications 
we have no prior information on p, and the same (lack of) prior information about each of the 
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components in the mixture. Thus, it makes sense to put a symmetric Dirichlet prior on the weights, 



and to take a prior on 6 that has the form 1X7=1 where tt : — > [0, 00) is a proper prior 

density on G. Let z = (z\,..., z m ) denote the observed data. It is well known that the resulting 



posterior density, tt(6, p\z), is intractable and highly multi-modal (see, for example, Jasra, Holmes 



and Stephens| 2005| ). Indeed, let E denote any one of the k\ permutation matrices of dimension k 



and note that w(0,p\z) = ir(E0, Ep\z). Thus, every local maximum of the posterior density has 
k\ — 1 exact replicas somewhere else in the parameter space. 

The standard DA algorithm for this mixture problem was introduced by Diebolt and Robert 



(1994) and is based on the following augmented model. Assume that {(Yi, Zi)Y^=\ we iid pairs 
such that Yi = j with probability pj, and, conditional on Y{ = j, Zi ~ Note that the 

marginal density of Zi under this two level hierarchy is just ([T]). Let y = (yi, . . . ,y m ) denote 
a realization of the YjS. The so-called complete data posterior density, ir((6,p), y\z), is just the 
posterior density that results when we combine our model for {(Yi, Zi)}™ =1 with the priors on p 
and 6 defined above. It is easy to see that 

5^7r((0,p),y|z) = n(e,p\z) , 
yeY 

where Y is the set of all sequences of length m consisting of integers from the set {1, . . . , k}. 
Hence, ir((0,p), y\z) can be used to build a DA algorithm as long as it is possible to sample from 
the conditionals, ir((6 ,p)\y, z) and ir(y\(0,p), z). We call it the mixture DA (MDA) algorithm. 
Note that the state space for the MDA chain is the Cartesian product of M fc ' and the fc-dimensional 
simplex, but |Y| = k m < 00. 

The MDA algorithm often converges very slowly because it moves between the symmetric 



modes of tt(6,p\z) too infrequently (Celeux, Hum and Robert 2000, Lee, Marin, Mengersen and 



Robert 2008| ). |Frtihwirth-Schnatter (2001 1 suggested adding a random label switching step to each 



iteration of the MDA algorithm in order to force movement between the modes. We show that 
the resulting Markov chain, which we call the FS chain, is a special case of the sandwich chain. 
Moreover, our theoretical results are applicable and imply that the spectrum of the operator defined 
by the FS chain dominates the spectrum of the MDA operator. To illustrate the extent to which the 
label switching step can speed convergence, we study two specific mixture models and compare the 
spectra associated with the FS and MDA chains. The first example is a toy problem in which we are 
able to get exact formulas for the eigenvalues. The second example is a normal mixture model that is 
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frequently used in practice, and we approximate the eigenvalues via classical Monte Carlo methods. 
The conclusions from the two examples are quite similar. Firstly, the MDA chain converges slowly 
and its rate of convergence deteriorates very rapidly as the sample size, m, increases. Secondly, 
the FS chain converges much faster and its rate does not seem as adversely affected by increasing 
sample size. 

The remainder of this paper is organized as follows. Section [2] is a brief review of the operator 
theory used for analyzing reversible Markov chains. Section [3] contains a string of results about the 
DA operator and its spectrum. Our main result comparing the DA and sandwich chains in the case 
where |Y| < oo appears in Section [4] Section [5] contains a detailed review of the MDA and FS 
algorithms, as well as a proof that the FS chain is a special case of the sandwich chain. Finally, 
in Section |6| the MDA and FS chains are compared in the context of two specific examples. The 
Appendix contains an eigen-analysis of a special 4x4 Mtm. 



2 Operator Theory for Reversible Markov Chains 

Consider the following generalized version of the problem described in the Introduction. Let X be a 
general space (equipped with a countably generated c-algebra) and suppose that fx : X — > [0, oo) is 
an intractable probability density with respect to the measure \i. Let p(x'\x) be a Mtd (with respect 
to n) such thatp(x'\x)fx(x) is symmetric in (x,x'), so the Markov chain defined by p is reversible 
with respect to fx(x). Assume that the chain is Harris ergodic, which means that it is irreducible, 
aperiodic and Harris recurrent ( Asmussen and Glynn||2010{ Meyn and Tweedie 1993 1. 



Define the Hilbert space 



L 2 {fx) = 1 5 : X -»• M : j^g 2 {x)f x {x)n{dx) < oo and j^g(x)f x (x)n(dx) = j , 



where inner product is defined as 

(9,h}= I g(x)h(x) fx(x) fi(dx) . 



x 



The corresponding norm is given by \\g\\ = (g, g). The Mtd p defines an operator P : L^fx) 
Lq(/x) that acts on g £ L^(fx) as follows: 

(Pg)(x) = / g{x)p{x'\x) /j,(dx') . 
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It is easy to show, using reversibility, that for g,h G Lq(/x), (Pg, h) = (g, Ph); that is, P is a 
self-adjoint operator. The spectrum of P is defined as 



MP) 



X G R : P - XI is not invertible 



}• 



There are two ways in which P — XI can fail to be invertible (Rudin 1991| Chapter 4). Firstly, 



P — XI may not be onto; that is, if there exists h G L^(fx) such that there is no g G L^(fx) for 
which ((P — XI) g) = h, then the range of P — XI is not all of Lq(/x), so P — XI is not invertible 
and A G Sp(P). Secondly, P — XI may not be one-to-one; that is, if there exist two different 
functions g,h G Lq(/x) such that ((P — XI)g) = ((P — XI)h), then P — XI is not one-to-one, so 
P-XIis not invertible and A G Sp(P). Note that, if ((P - XI) g) = ((P- Ai», then Pg* = Xg* 
with g* = g — h, and A is called an eigenvalue with eigen-function g*. We call the pair (A, g*) an 
eigen-solution. 

Let Lq i(fx) denote the subset of functions in Lq(/x) that satisfy f x g 2 (x) fx{x) /J,(dx) = 1. 
The (operator) norm of P is defined as 

||P|| = sup \\Pg\\ . 

A simple application of Jensen's inequality shows that the nonnegative quantity ||P|| is bounded 
above by 1. The norm of P is a good univariate summary of Sp(P). Indeed, define 



l P = inf {Pg, g) 



and up = sup (Pg,g) . 



It follows from standard linear operator theory that inf Sp(P) = lp, sup Sp(P) = up, and ||P|| = 
max {— lp, up}. Consequently, Sp(P) C [— ||P||, ||P||] C [—1, 1]. Another name for ||P|| in this 
context is the spectral radius, which makes sense since ||P|| represents the maximum distance that 
Sp(P) extends away from the origin. The quantity 1 — ||P|| is called the spectral gap. 

It is well known that ||P|| is closely related to the convergence properties of the Markov chain 
defined by p ( |Liu, W ong an d Kong| 1995[ Rosenthal 2003 1. In particular, the chain is geometrically 



ergodic if and only if ||P|| < 1 (Roberts and Rosenthal 1997). There is an important practical 



advantage to using an MCMC algorithm that is driven by a geometrically ergodic Markov chain. 
Indeed, when the chain is geometric, sample averages satisfy central limit theorems, and these allow 



for the computation of asymptotically valid standard errors for MCMC-based estimates (Flegal, 
Hara n and Jones| |2008t |Jones, Haran, Caffo and Neath] |2006[ ). We note that geometric ergodicity 
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of reversible Monte Carlo Markov chains is typically not proven by showing that the operator norm 



is strictly less than 1, but rather by establishing a so-called geometric drift condition (Jones and 
Hobert[|200T 



If |X| < oo, then P is simply the Mtm whose (i,j)th element is p(j\i), the probability that the 



chain moves from i to j. In this case, Sp(P) is just the set of eigenvalues of P (see e.g. Mira and 
GeyerJ 1999 1. The reader is probably used to thinking of 1 as an eigenvalue for P because P satisfies 
the equation PI = 1, where 1 denotes a vector of ones. However, the only constant function 
in Lg is the zero function, so (1,1) is not a viable eigen-solution in our context. Furthermore, 
irreducibility implies that the only vectors v that solve the equation Pv = v are constant. It follows 
that 1 ^ Sp(P). Aperiodicity implies that —1 ^ Sp(P). Hence, when X is a finite set, ||P|| is 
necessarily less than one. In the next section, we return to the DA algorithm. 

3 The Spectrum of the DA Chain 

Suppose that Y is a second general space and that u is a measure on Y. Let / : X x Y — > [0, oo) be 
a joint probability density with respect to \i x v. Assume that f Y f(x, y) v(dy) = fx(x) and that 
simulating from the associated conditional densities, fx\Y{'\y) an d fY\x('\ x )> * s straightforward. 
(For convenience, we assume that fx and fy are strictly positive on X and Y, respectively.) The 
DA chain, {X„}^L , has Mtd (with respect to f£) given by 

k(x'\x) = J f x \Y(x'\y) fY\x(y\ x ) v(dy) . (2) 

It is easy to see that k{x'\x) fx{x) is symmetric in (x, x'), so the DA chain is reversible with respect 
to fx- We assume throughout this section and the next that all DA chains (and their conjugates) 
are Harris ergodic. (See Hobert| ( |2011 1 for a simple sufficient condition for Harris ergodicity of the 



DA chain.) If the integral in ((2]) is intractable, as is nearly always the case in practice, then direct 
simulation from k(-\x) will be problematic. This is why the indirect two-step procedure is used. 
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Liu et al. (1994) showed that the DA chain satisfies an important property that results in a 
positive spectrum. Let K denote the operator defined by the DA chain. For g 6 L^(fx), we have 



(Kg,g) 



(Kg)(x) g(x) f x (x) p,(dx) 



g(x') k(x'\x) n(dx') 



g{x) fx(x) (Jt(dx) 



fx\y{x'\y) fy\x{y\x) v{dy) 



fi(dx') 



g(x) fx(x) fi(dx) 



g{x) fx\v{x\y) fi(dx) 



f Y {y)v(dy) > 0, 



which shows that K is a positive operator. It follows that Ik > 0, so Sp(K) C [0, \\K\\] C [0, 1] 
and ||i^|| = supSp(i^). 

In most applications of the DA algorithm, fx is a probability density function (with respect to 
Lebesgue measure), which means that X is not finite. Typically, when |X| = oo, it is difficult to 
get a handle on Sp(K), which can be quite complex and may contain an uncountable number of 
points. However, if K is a compact operatoiQ then Sp(iC) has a particularly simple form. Indeed, 
if | X | = oo and K is compact, then the following all hold: (i) the number of points in Sp(Zf) is at 
most countably infinite, (ii) {0} G Sp(Zf), (iii) {0} is the only possible accumulation point, and (iv) 
any point in Sp(K) other than {0} is an eigenvalue. In the remainder of this section we prove that, 
if | X | = oo and |Y| = d < oo, then K is a compact operator and Sp(i^) consists of the point {0} 
along with d — 1 eigenvalues, and these are exactly the d — 1 eigenvalues of the Mtm that defines 
the conjugate chain. It follows immediately that the DA chain is geometrically (in fact, uniformly) 
ergodic. Moreover, K has a finite spectral decomposition that provides very precise information 



about the convergence of the DA chain (Diaconis et al. 2008). Indeed, let {(A^, <7i)}f =1 denote a 
set of (orthonormal) eigen-solutions for K. If the chain is started at Xq = x, then the x 2 -distance 



between the distribution of X n and the stationary distribution can be expressed as 

|2 d-1 



\k n (x'\x) - f x (x')\ 
fx(x') 



n{dx') = Y J ^ n gKx) 



(3) 



i=l 



where k n {-\x) is the n-step Mtd; that is, the density of X n given Xq = x. Of course, the x 2 -distance 



is an upper bound on the total variation distance (see, for example, |Liu et al.] [1995 ). Since the A,s 



'The operator K is defined to be compact if for any sequence of functions Qi in L D (/x) with \\gi\\ < 1, there is a 
subsequence gi . such that the sequence Kgi j converges to a limit in Lq(/x). 
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are the eigenvalues of the Mtm of the conjugate chain, there is some hope of calculating, or at least 
bounding them. 

Let Lq(/y) be the set of mean-zero, square integrable functions with respect to fy. In a slight 
abuse of notation, we will let (•, •) and ||-|| do double duty as inner product and norm on both 
Lq(/x) and on Lq(/y). We now describe a representation of the operator K that was developed 



and exploited by Diaconis et al. (2008l (see also Buja 1990). Define Q : L^(fx) — > L^Uy) and 
Q* '■ l oUy) -> Llifx) as follows: 

(Qg)(y)= / 9(x) f x \Y(x\y) fi(dx) and (Q*h)(x)= / h(y) f Y \ x (y\x)v{dy) . 
Jx Jy 



Note that 



(Qg,h)= / {Qg){y)h{y)fy{y)v{dy) 



g{x) fx\Y(x\y) n{dx) 



x 



x 

g( x ) 



h(y) fY{y)v(dy) 



h(v) fY\x(y\x) v{dy) 



fx(x)n{dx) 



(g,Q*h) , 



which shows that Q* is the adjoint of Q. (Note that we are using the term "adjoint" in a somewhat 
non-standard way since (Qg, h) is an inner product on L^fy), while (g, Q*h) is an inner product 
on Lq(/x).) Moreover, 



(Kg)(x)= / g(x') k(x'\x) fi(dx') 



g( x ') 



fx\Y(x'\y) fy\x(y\x) v(dy) 



fj,(dx) 



g{x') fx\Y{x'\y) n(dx') 



{Qg){y) fY\x{y\x)v{dy) 



((Q*Q)g)(x) , 



fv\x(y\x) v{dy) 



which shows that K = Q*Q. As in Section [TJ consider the conjugate Markov chain whose Mtd 
(with respect to u) is given by 

Hy'\y) = / fY\x(y'\x) f X \Y(x\y) Kdx) . (4) 

Jx 
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Obviously, k(y'\y) is reversible with respect to fy. Furthermore, it is easy to see that K = QQ* , 
where K : Lg(/y) — > L^(fy) is the operator associated with k. 

Now suppose that (A, g) is an eigen-solution for K; that is, (Kg)(x) = Xg(x), or, equiva- 
lent^, ((Q*Q)g)(x) = Xg(x). Applying the operator Q to both sides yields, {Q((Q*Q)g)){y) = 
X(Qg)(y), but we can rewrite this as (K(Qg))(y) = X(Qg)(y), which shows that (X,Qg) is an 



eigen-solution for K. (See Buja (1990) for a similar development.) Of course, the same argument 
can be used to convert an eigen-solution for K into an eigen-solution for K. We conclude that K 
and K share the same eigenvalues. Here is a precise statement. 

Proposition 1. If (X,g) is an eigen-solution for K, then (X,(Qg)) is an eigen-solution for K. 
Conversely, if(X, h) is an eigen-solution for K, then (A, (Q*h)) is an eigen-solution for K. 



Remark 1. Diaconis et al. \2008 \ describe several examples where the eigen-solutions of K and 
K can be calculated explicitly. These authors studied the case where fx\Y( x \y) a univariate 
exponential family (with y playing the role of the parameter), and fy{y) is the conjugate prior. 



The next result, which is easily established using minor extensions of results in Retherford 



( 1993 ) Chapter VII, shows that compactness is a solidarity property for K and K. 



Proposition 2. K is compact if and only if K is compact. 

Here is the main result of this section, which relates the spectrum of the DA chain to the spec- 
trum of the conjugate chain. 

Proposition 3. Assume that |X| = oo and |Y| = d < oo. Then K is a compact operator and 
Sp(K) = {0}USp(K). 

Proof. Since |Y| < oo, K is a compact operator. It follows from Proposition [2] that K is also 
compact. Hence, {0} G Sp(K), and aside from {0}, all the elements of Sp(K) are eigenvalues of 
K. But we know from Propositionlllthat K and K share the same eigenvalues. □ 



Remark 2. 



Liu et al. 



's (1994) Theorem 3.2 states that \\K\\ = \\K\\ (regardless of the cardinalities 
ofX andY). Proposition^can be viewed as a refinement of this result in the case where | Y| < 00. 
See also \Roberts a nd Rosenthal ( 2001 ). 



In the next section, we use Proposition [3] to prove that the spectrum of the sandwich chain 
dominates the spectrum of the DA chain. 
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4 Improving the DA algorithm 

Suppose that R(y, dy') is a Markov transition kernel on Y that is reversible with respect to fy (y). 
Let {Xn\^ = Q be the sandwich chain on X whose Mtd is given by 



k{x'\x) = J j f xlY (x'\y') R(y, dy') f Y]x {y\x) u{dy) . (5) 

Again, routine calculations show that the sandwich chain remains reversible with respect to the 
target density fx- Moreover, if we can draw from R(y, •), then we can draw from k(-\x) in three 
steps. First, draw Y ~ fy\x{'\ x )> ca U the result y, then draw Y' ~ R(y, •), call the result y', and 
finally draw X' ~ f x \y i'W)- 

Note that k is not defined as the integral of the product of two conditional densities, as in Q. 
However, as we now explain, if R satisfies a certain property, called idempotence, then k can be 
re-expressed as the Mtd of a DA chain. The transition kernel R(y, dy') is called idempotent if 
R 2 (y, dy') = R(y, dy') where R 2 (y, dy') = f y R(y, dw) R(w, dy'). This property implies that, if 
we start the Markov chain (defined by R) at a fixed point y, then the distribution of the chain after 
one step is the same as the distribution after two steps. For example, if R(y, dy') does not depend 
on y, which implies that the Markov chain is just an iid sequence, then R is idempotent. Here is a 
more interesting example. Take Y = R and R(y, dy') = r(y'\y) dy' with 

r(y'\y) = e" l?/l i[o,oo)(y)i[o,oo){y') + i(-oo,o){y)i(-oo,o){y') 

It is easy to show that f R r(y'\w) r(w\y) dw = r(y'\y), so R is indeed idempotent. Note that the 
chain is reducible since, for example, if it is started on the positive half-line, it can never get to the 
negative half-line. In fact, reducibility is a common feature of idempotent chains. Fortunately, the 
sandwich chain does not inherit this property. 

Hobert and Marchev] p008) proved that if R is idempotent, then 



k(x'\x) = I f xlY (x'\y) f Ylx (y\x) u{dy) , (6) 

where 



Y 



f*{x,y) = f Y (y) / fx\y(x\y') R(y,dy') 



Y 



Note that /* is a probability density (with respect to fj, X v) whose x and y-marginals are fx 
and fy- What is important here is not the particular form of /*, but the fact that such a density 
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exists, because this shows that the sandwich chain is actually a DA chain based on the joint density 
f*(x, y). Therefore, we can use the theory developed in Section [5] to analyze the sandwich chain. 
Let K : Lq(Jx) — > L"l{fx) denote the operator defined by the Mtd k. 
Corollary 1 states that \\K\\ < \\K\\ (see also 



Hobert and Marchev 



s (2008 l 



Hobert and Roman 



2011). Here is a refinement of 



that result in the case where | Y| < oo. 

Theorem 1. Assume that |X| = oo, |Y| = d < oo and that R is idempotent. Then K and K are 
both compact operators and each has a spectrum that consists exactly of the point {0} and d — 1 
eigenvalues in [0, 1). Furthermore, if we denote the eigenvalues of K by 



< X d -i < \ d -2 < ■ • • < Ai < 1 , 



and those of K by 



< X d -i < X d - 2 < • • • < Ai < 1 , 
then Aj < Xi for each i G {1, 2, . . . ,d— 1}. 

Proof Since R is idempotent, both k and k are DA Markov chains. Moreover, in both cases, the 
conjugate chain lives on the finite space Y, which has d elements. Therefore, Proposition [3] implies 
that K and K are both compact and each has a spectrum consisting of the point {0} and d — 1 
eigenvalues in [0, 1). Now, Corollary 1 of 



Hobert and Marchev 



(20081 implies that K - K is a 



positive operator. Thus, for any g G Lq(/x), 



(Kg,g) < (Kg,g) 



(9,9) 



(9,9) 



The eigenvalue ordering now follows from an extension of the argument used to prove Mira and 
Geyer 's ( 1999) Theorem 3.3. Indeed, the Courant-Fischer-Weyl minmax characterization of eigen- 
values of compact, self-adjoint operators (see, e.g., Voss| 2003 ) yields 

(Kg, 9) / . (Kg, 9) 



Xi = min max — U1,J ' < mm max 

dim(V)=i-i g&v 1 - ,<?^o {g,g) dim(V)=i-i gev ± , g ^o {g,g) 



Xi , 



where V denotes a subspace of L^(fx) with dimension dim(y), and V 1 - is its orthogonal comple- 



ment. 



□ 



Theorem[l]shows that, unless the two spectra are exactly the same, Sp(i^T) is closer than Sp(ET) 
to the ideal spectrum, {0}. In fact, in all of the numerical comparisons that we have performed, it 
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has always turned out that there is strict inequality between the eigenvalues (except, of course, when 
they are both zero). When the domination is strict, there exists a positive integer N such that, for all 

n > N, 

f \~k n (x'\x)-f x (x')\ 2 f \k n (x'\x)-f x (x')\ 2 

/ TTT\ ^ dx < / TTT\ V(dx ) . 

Jx fx{x') Jx Jx(x') 

Indeed, let {(Xi, gi)}fli denote a set of (orthonormal) eigen-solutions of K. Then, according to 

([3]), the x 2 -distance between the distribution of X n and the stationary distribution is given by 

d-l 

E^ 2 (x). (7) 

i=l 

Now, fix i G {1, . . . , d — 1}. If X{ = Xi = 0, then the ith term in the sum is irrelevant. On the other 
hand, if < X\ < Xi, then, no matter what the values of gi{x) and gi(x) are, Xj n gf(x) will be less 
than X? n gf(x) for all n eventually. 

In the next section, we provide examples where the sandwich chain converges much faster than 
the DA chain, despite the fact that the two are essentially equivalent in terms of computer time per 
iteration. 

5 Improving the DA Algorithm for Bayesian Mixtures 
5.1 The model and the MDA algorithm 

Let C MJ and consider a parametric family of densities (with respect to Lebesgue or counting 
measure on W) given by {hg(-) : 6 £ 0}. We work with a ^-component mixture of these densities 
that takes the form 

k 

f(z\e, P ) = Y,pM(z) , (8) 

where = (Q\, . . . ,0k) T £ © fe andp = (pi, . .. ,Pk) T 6 Sfc, where 

S k := {pGl fc : K 6 [0, 1] and Pl + ■ ■ ■ +p k = l} . 

Let Z\, . . . , Z m be a random sample from / and consider a Bayesian analysis of these data. We 
take the prior for 9 to be J^Li n(9j), where vr : 6 -> [0, oo) is a proper prior density on ©. The 
prior on p is taken to be the uniform distribution on Sk- (The results in this section all go through 
with obvious minor changes if the prior on p is taken to be symmetric Dirichlet, or if p is known 
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and all of its components are equal to 1/k.) Letting z = (z\, . . . , z m ) denote the observed data, the 
posterior density is given by 

(k-i)iis k (p)\nU< e j)]f( z \ e >p) 



n(0,p\z) 



m(z) 
where 

m 

f(z\0,p)=H 



(9) 



i=l 



k 

i=i 



and m(z) denotes the marginal density. The complexity of this posterior density obviously depends 
on many factors, including the choices of hg and ir, and the observed data. However, the versions of 
tt(0,p\z) that arise in practice are nearly always highly intractable. Moreover, as we now explain, 
every version of this posterior density satisfies an interesting symmetry property, which can render 
MCMC algorithms ineffectual. 

The prior distribution on (6, p) is exchangeable in the sense that, if E is any permutation matrix 
of dimension k, then the prior density of the point (0, p) is equal to that of (EG, Ep). Furthermore, 
the likelihood function satisfies a similar invariance. Indeed, f(z\EO, Ep) does not vary with E. 
Consequently, ir(EO, Ep\z) is invariant to E, which means that any posterior mode has k\ — 1 exact 
replicas somewhere else in the space. Now, if a set of symmetric modes are separated by areas of 
very low (posterior) probability, then it may take a very long time for a Markov chain (with invariant 
density n(0, p\z)) to move from one to the other. 

We now describe the MDA algorithm for exploring the mixture posterior. Despite the fact that 



this algorithm has been around for many years (Diebolt and Robert 1994| ), we provide a careful 
description here as this will facilitate our development of the FS algorithm. Consider a new (joint) 
density given by 

k 

f(z,y\e, P ) = Y,Pjh3}(y) h oA z ) ■ ( 10 ) 

Integrating z out yields the marginal mass function of Y, which is Ylj=i Pjl{j}(v)- Hence, Y is a 
multinomial random variable that takes the values 1, . . . , k with probabilities pi, . . . ,Pk- Summing 
out the y component leads to 

k k 

^2f(z,y\e,p) = Y t p j h .(z), (11) 
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which is just ([8]>. Equation (111 establishes Y as a latent variable. Now suppose that {(Yi, 
are iid pairs from ( [T0| ). Their joint density is given by 

k 



f(z,y\0,p) = Y\ 



i=l 



3=1 



where y = (yi, . . . , y m ) takes values in Y, the set of sequences of length m consisting of positive 
integers between 1 and k. Combining f(z,y\0,p) with our prior on (9,p) yields the so-called 
complete data posterior density given by 

(k-i)\i Sh (p)[nU< e i)]f( z >v\ >p) 

m{z) 

This is a valid density since, by ( [TT] ), 



n(0,p,y\z) 



(12) 



which in turn implies that 



^2f(z,y\e,p) = f(z\e, P ) 



^2ir(9,p,y\z) = 7r{6,p\z) . 



(13) 



In fact, ( pL3[ > is the key property of the complete data posterior density. In words, when the y 
coordinate is summed out of n(6,p, y\z), we are left with the target density. Hence, we will have 
a viable MDA algorithm as long as straightforward sampling from n(0, p\y, z) and n(y\6, p, z) is 
possible. Note that the roles of x and y from Sections[T||3]and|4]are being played here by (6, p) and 
y, respectively. 



Now consider sampling from the two conditionals. First, it follows from ( 12 1 that 



TT(y\6,p,z) = YI 



i=l 



Ylj=iVjI{j}{yi)he 3 {zi) 



(14) 



Therefore, conditional on (0,p, z), the YjS are independent multinomial random variables and Yi 
takes the value j with probability Pjhe^Zi)/ ( Ef=iM^fi( z i)) f° r J e {!>•••> Consequently, 
simulating from n(y\6, p, z) is simple. 

A two-step method is used to sample from n(0,p\y, z). Indeed, we draw from n(p\y, z) and 
then from ir(0\p, y, z). It follows from ( fT2] ) that 

k 

7r(p\0, y, z) oc I Sk (p) J[pj j , 

3=1 
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where cj = I{j}(Vi)- This formula reveals two facts: (i) given (z,y), p is conditionally 

independent of 6, and (ii) the conditional distribution of p given (z, y) is Dirichlet. Thus, it is easy 
to draw from n(p\y, z), and our sequential strategy will be viable as long as we can draw from 
ir(6\p, y, z). Our ability to sample from n(0\p, y, z) will depend on the particular forms of hg and 
the prior -k. In cases where it is a conjugate prior for the family fig, it is usually straightforward to 



draw from n(G\p, y, z). For several detailed examples, see Chapter 9 of Robert and Casella (2004l. 



The state space of the MDA chain is X = fc x and its Mtd is given by 

k(e',p'\e, P ) = Y,Ao',p'\y,z)TT(y\e, P ,z) . 

yev 

Since |Y| = k m , Proposition [5] implies that the operator K : Lq[-k{6 ,p\z)) — > L\ (vr(0, p\z)) 
defined by k(6',p'\6, p) is compact and 

Sp(-KT) = {0, Afcm_ l5 A / fcm_ 2 , . . . , Al} , 

where < Afcm_i < Afcm_ 2 < • • • < Ai < 1, and the A»s are the eigenvalues of the k m x k rn Mtm 
defined by 

Hy'\y)= / ■n(y'\0,p,z)ir(6,p\y,z)dpd0 . 

Je k Js k 

As far as we know, there are no theoretical results available concerning the magnitude of the AjS. 
On the other hand, as mentioned in Section[T| there is a great deal of empirical evidence suggesting 
that the MDA chain convergences very slowly because it moves between the symmetric modes of 
the posterior too infrequently. In the next section, we describe an alternative chain that moves easily 
among the modes. 



5.2 Fruhwirth-Schnatter's algorithm 

One iteration of the MDA chain can be represented graphically as (9,p) — >• y — > (9',p'). To 
encourage transitions between the symmetric modes of the posterior, Fruhwirth-Schnatter (2001) 
suggested adding an extra step to get (6, p) — » y — ^ y' — » (0' ,p'), where the transition y — > y' is a 
random label switching move that proceeds as follows. Randomly choose one of the k\ permutations 
of the integers 1, . . . , k, and then switch the labels in y according to the chosen permutation to get 
y' . For example, suppose that m = 8, k = 4, y = (3,3,4,1,3,3,4,3), and that the chosen 
permutation is (1324). Then we move from y to y' = (2, 2, 1, 3, 2, 2, 1, 2). Using both theory and 
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examples, we will demonstrate that Friihwirth-Schnatter s (2001 ) Markov chain, which we call the 
FS chain, explores ir(0,p\z) much more effectively than the MDA chain. 

To establish that the results developed in Section |4] can be used to compare the FS and MDA 
chains, we must show that the FS chain is a sandwich chain with an idempotent r. That is, we must 
demonstrate that the Mtd of the FS chain can be expressed in the form 

k(0',p'\G,p) = J2J2 K(e',p'\y',z)r(y'\y)ir(y\e,p,z) , (15) 

where r(y'\y) is a Mtm (on Y) that is both reversible with respect to 

n{y\z)= / Tr(0,p,y\z)dOdp , 

and idempotent. We begin by developing a formula for r(y'\y). Let &k denote the set (group) of 
permutations of the integers 1, . . . , k. For a G (3&, let ay represent the permuted version of y. 
For example, if y = (3, 3, 4, 1, 3, 3, 4, 3) and a = (1324), then ay = (2, 2, 1, 3, 2, 2, 1, 2). The 
label switching move, y — > y', in the FS algorithm can now be represented as follows. Choose a 
uniformly at random from &k and move from y to y' = ay. Define the orbit of y G Y as 

Oy = {y' ^ Y ■ y' = a y f° r some a G &f.\ . 

The set O y simply contains all the points in Y that represent a particular clustering (or partitioning) 
of the m observations. For example, the point y = (3, 3, 4, 1, 3, 3, 4, 3) represents the clustering of 
the m = 8 observations into the three sets: {1, 2, 5, 6, 8}, {3, 7}, {4}. And, for any a G &k, &y 
represents that same clustering because all we're doing is changing the labels. 

We now show that, if y is fixed and a is chosen uniformly at random from <£>&, then the random 
element ay has a uniform distribution on O y . Indeed, suppose that y contains u distinct elements, 
so u G {1, 2, . . . , k}. Then, for any fixed y' G O y , exactly (k — u)\ of the k\ elements in satisfy 
ay = y'. Thus, the probability that ay equals y' is given by (k — u)\/k\, which does not depend 
on y'. Hence, the distribution is uniform. (Note that this argument implies that \O y \ = k\/[k — u)\, 
which can also be shown directly.) Therefore, we can write the Mtm r as follows: 

r(y'\y) = jOyihoyyiy') ■ 

Since the chain driven by r cannot escape from the orbit (clustering) in which it is started, it is 
reducible. (Recall from Section|4]that reducibility is a common characteristic of idempotent Markov 
chains.) 
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A key observation that will allow us to establish the reversibility of r is that n(y\z) = ir(ay\z) 
for all y G Y and all a G &k- Indeed, 



3=1 



dp )dG . 



Let ay = y' = (y[, . . . , y' m ). Now, since y- = a(j) y { = j, we have 

k k 
3=1 3=1 



Hence, 



ir(ay\z) 



(k-l)l 
m(z) 



m r k 



y J s k i= i l 7=1 j j 



The fact that ir(y\z) = n{ay\z) can now be established through a couple of simple arguments 
based on symmetry. 

We now demonstrate that the Mtm r satisfies detailed balance with respect to Tr(y\z); that is, 
we will show that, for any y,y' G Y, r(y'\y) ir(y\z) = r(y\y')ir(y'\z). First, a little thought 
reveals that, for any two elements y and y', only one of two things can happen: either O v = O y i or 

Oy fl Oy, = 0. If Oy D O >y, = ^, toOl J {0 „}(y') = I {O y ,}{V) = ^ ^ T \v) = T^rf) = Mid 

detailed balance is satisfied. On the other hand, if O y = O y i, then I{o y }{y') = I{o ,}{y) = 1 an d 
Vl^yl = Vl^y'l' so r {y'\y) = r (y\y')> an d tne common value is strictly positive. But y' G O y 
implies that y' = ay for some a G <3&. Thus, ir{y\z) = iv(y'\z), and detailed balance holds. 

Finally, it is intuitively clear that r is idempotent since, if we start the chain at y, then one step 
results in a uniformly chosen point from O y . Obviously, the state after two steps is still uniformly 
distributed over O y . Here's a formal proof that r 2 {y'\y) = r{y'\y). For y, y' G Y, we have 



\y'\y) = ^2 r(y'\w)r(w\y) 

_,, \ ( - y iii\ \ K - / v\ 



= To~\ Yl w-\ho w }(y) 

|U ™ weo y lUwl 



\°„, 

r(y'\y) , 



weOy 
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where the fourth equality follows from the fact that w £ O y => O w = O y . 

We have now shown that the Mtd of the FS chain can indeed be written in the form ( [13] ) with 
an appropriate r that is reversible and idempotent. Hence, Theorem[T]is applicable and implies that 
the operators defined by the two chains are both compact and each has a spectrum consisting of the 
point {0} and k m — 1 eigenvalues in [0, 1). Moreover, Aj < Aj for each i € {1,2,..., k m — 1}, 
where {Ai}^ -1 and {Aj}k =1 -1 denote the ordered eigenvalues associated with the FS and MDA 
chains, respectively. 

Interestingly, in the special case where m = 1, the FS algorithm actually produces an iid se- 
quence from the target distribution. Recall that ir(y\z) = n(ay\z) for all y G Y and all a £ &k- 
Thus, all the points in O y share the same value of tt(-\z). When m = 1, Y contains only k points 
and they all exist in the same orbit. Thus, ir(y\z) = 1/k for all y £ Y. Moreover, since there is 
only one orbit, r(y'\y) = 1/k for all y' £ Y; i.e., the Markov chain corresponding to r is just an 
iid sequence from the uniform distribution on Y. In other words, the label switching move results 
in an exact draw from ir(y'\z). Now recall the graphical representation of one iteration of the FS 
algorithm: (6,p) — » y — ^ y' — » (6',p'). When m = 1, the arguments above imply that, given 
(0,p), the density of (y, y', 0',p') is 

n(y\0,p, z)r(y'\y)ir(0' ,p'\y' , z) = ir(y\6,p, z)ir(y'\z)ir(0' ,p'\y' , z) . 

Thus, conditional on (6, p), y and (y 1 , 6 1 , p') are independent, and the latter has density 

7r(y'\z)iT(0',p'\y',z) = k(6' ,p' ,y'\z) . 

It follows that, marginally, (6',p') ~ tt(6' ,p'\z), so the FS algorithm produces an iid sequence 
from the target posterior density. When m = 1, |Y| = k m = k. Thus, while the spectrum of the 
MDA operator contains k — 1 eigenvalues, at least one of which is strictly positive, the spectrum of 
the FS operator is the ideal spectrum, {0}. 

In the next section, we consider two specific mixture models and, for each one, we compare the 
spectra associated with FS and MDA chains. The first example is a toy problem where we are able 
to get exact formulas for the eigenvalues. The second example is a normal mixture model that is 
frequently used in practice, and we approximate the eigenvalues via classical Monte Carlo methods. 
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6 Examples 



6.1 A toy Bernoulli mixture 

Take the parametric family hg to be the family of Bernoulli mass functions, and consider a two- 
component version of the mixture with known weights both equal to 1/2. This mixture density 
takes the form 

f(z\r, S ) = ^r z (l-r) 1 -* + ±s*(l-s) 1 -* , 

where z G {0, 1} and = (r, s). To simplify things ever further, assume that r, s G {p, 1 — p} 
where p G (0, 1/2) is fixed; that is, the two success probabilities, r and s, can only take the values 
p and 1 — p. Hence, (r, s) G X = {(p, p), (p, 1 — p), (1 — p, p), (1 — p, 1 — p)}. Our prior for 
(r, s) puts mass 1/4 on each of these four points. A simple calculation shows that the posterior mass 
function takes the form 

ir(r,s\z) = — — —. — — , 

2 m p m i(l - p) m ~ m i + 2 m p m - m i(l - p) m i + 2 

where z = (z\, . . . , z m ) G {0,l} m denotes the observed data, and m\ denotes the number of 
successes among the m Bernoulli trials; that is, m\ = XX=i z, i- While we would never actually use 
MCMC to explore this simple four-point posterior, it is both interesting and useful to compare the 
FS and MDA algorithms in this context. 

As described in Subsection |5.1| the MDA algorithm is based on the complete data posterior 
density, which is denoted here by ir(r, s, y\z). (The fact that p is known in this case doesn't really 
change anything.) Of course, all we really need are the specific forms of the conditional mass 



functions, 7r(y\r, s, z) and 7r(r, s\y, z). It follows from the general development in Subsection 5.1 



that, given (r,s,z), the components of y = (yi, y2, ■ ■ ■ , y m ) are independent multinomials with 
mass functions given by 

I {l} { yi )r^{l-r) l -^+I {2] ( yi )s^(l- S ) 1 -^ 
ir{yi\r, s, z) = — — r— — r— . 

7* * ( 1 — V) ~\~ S^M 1 — Sj ^ 

Furthermore, it is easy to show that, given (y,z), r and s are independent so 7r(r, s\y,z) = 
7r(r\y, z) ir(s\y, z). Now, for j G {1, 2} and k G {0, 1}, let rrijk denote the number of (yj, Zi) 
pairs that take the value (j, k). (Note that min + mn = c\ and raw + ^21 = mi.) Then we have 

/ {p} (r)p m "(l - p) m ^ + / {1 _ p} (r)p"^(l - p) m " 



7r(r\y,z) 



p"Hi(l — p\m w -\- p m w{\ — p) m n 
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and 



w(s\y,z) 



I {p} (s)p^ (1 - p) m ™ + I {1 _ p} ( S ) P ™™(l - p)™ 



p™2i(\ — p^m 20 _|_ p™2o(l _ p)"i 2 i 

The state space of the MDA chain is X = { (p, p), (p, 1 — p), (1 — p, p), (1 — p, 1 — p)}, which has 
only four points. Hence, in this toy Bernoulli example, we can analyze the MDA chain directly. Its 
Mtm is 4 x 4 and the transition probabilities are given by 



k(r', s'\r, s) = 7r ( r '' s '|j/> z ) n(v\ r i s -> z ) 
veY 



(16) 



where Y = {1, 2} m . We now perform an eigen-analysis of this Mtm. Note that 7r(r', s'|y, z) and 
7r(y|r, s, z) depend on y only through m\o, m,\\, 777,20 and 77721. If we let = m — m\, then we 
can express the transition probabilities as follows: 

7 w (rV (1 " P) 3 + / { i_ p} (rV(l - p) 4 



fc(/,s'|r,s) = 

j=0 j=0 



mi 



m 
J 



pi(l_p)j+pj(l_p)i 



pmi-i^ _ p)m -j _|_ p»"0-i(l — p) m l-» 

Now, for /c = 0, 1, 2 define 



mi- 



j „rai-i 



m -J 



(r + s) m i(2 - r - s) 



m 



W k {p) = 
mi mo 

EE 

i=0 j=0 



mi 

i 



mo 
3 



pk(m ~j+i) _ p\k(m\-i+j) 



(p*(l - p)J + p?'(l - p)*) (p™!-^! - p)™0-j + pmo-j(l - p)" 1 !-*) 



Using this notation, we can write the Mtm as follows: 

p m l(l-p) m /x 1 / n 1 / \ p m 0(l-p) m l , \ 

pmi(l _ p) m ° Wl (p) W 2 { P ) P m (l - P ) m Wo(p) p m °(l ~ p) mi Wi{ P ) 

pmi(l _ p)^0 Wl {p) p m {l- p) m Wo{p) W 2 {p) p m °(l ~ p) mi Wi{p) 

Wo{p) ^nWi(p) ^awi(p) p ° ( ^ p) - wo(p) 



k 



l(l-p) m 0_ 



We have ordered the points in the state space as follows: (p, p), (p, 1— p), (1— p, p), and (1— p, 1— p). 
So, for example, the element in the second row, third column is the probability of moving from 
(p, 1 — p) to (1 — p, p). Note that all of the transition probabilities are strictly positive, which 
implies that the MDA chain is Harris ergodic. 

Of course, since k is a Mtm, it satisfies kvo = Xqvo where vq = 1 and Ao = 1. Again, (vo, Ao) 
does not count as an eigen-solution for us because we are using L^(fx) instead of L 2 (fx), and the 
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only constant function in Lq(/x) is 0. For us, there are three eigen-solutions, and we write them as 
(v{, Aj), i G {1, 2, 3}, where < A3 < A2 < Ai < 1. Note that the first and fourth rows of k are 
identical, which means that A3 = 0. The remaining eigen-solutions follow from the general results 
in the Appendix. Indeed, 

\i = Mp)-p m 0--p) m Mp), 

and the corresponding eigen- vector is v\ = (0, 1, —1, 0) T . Finally, 

g(p)wo(p) 



A 2 



g(p)wi(p) 



and v 2 = (a, 1, 1, a) T , where g(p) = p mi (l - p) m ° + p m °(l - p) mi and 

g(p)w (p) - 2 m 



a = 



2 m g(p)w 1 {p) 

(The fact that A2 < Ai actually follows from our analysis of the FS chain below.) We now use these 
results to demonstrate that the MDA algorithm can perform quite poorly for the Bernoulli model. 

Consider a numerical example in which m = 10, p = 1/10 and the data are z\ = • • • = z§ = 
and zq = ■ ■ ■ = zio = 1. The posterior mass function is as follows: 

ir(p,p\z) = 7r(l-p, l-p\z) = 0.003 and n(p, 1 - p\z) = tt(1 - p, p\z) = 0.497 . 

So there are two points with exactly the same very high probability, and two points with exactly the 
same very low probability. The MDA chain converges slowly due to its inability to move between 
the two high probability points. Indeed, the Markov transition matrix in this case is: 

0.10138 0.39862 0.39862 0.10138 

0.00241 0.99457 0.00061 0.00241 

0.00241 0.00061 0.99457 0.00241 

0.10138 0.39862 0.39862 0.10138 

Suppose we start the chain in the state (p, 1 — p). The expected number of steps before it reaches 
the other high probability state, (1 — p, p), is quite large. First, we expect the chain to remain in the 
state (p, 1 — p) for about 1/(1 — 0.99457) w 184 iterations. Then, conditional on the chain leaving 
(p, 1 — p), the probability that it moves to (p, p) or (1 — p, 1 — p) is about 0.89. And if it does 
reach (p, p) or (1 — p, 1 — p), there is still about a 40% chance that it will jump right back to the 
point (p, 1 — p), where it will stay for (approximately) another 184 iterations. All of this translates 
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into slow convergence. In fact, the two non-zero eigenvalues are (Ai, A2) = (0.99395, 0.19795). 
Moreover, the problem gets worse as the sample size increases. For example, if we increase the 
sample size to m = 20 (and maintain the 50:50 split of 0s and Is in the data), then (Ai, A2) = 
(0.99996, 0.15195). Figure[T]shows how the dominant eigenvalue, Ai, changes with sample size for 
several different values of p. We conclude that, for fixed p, the convergence rate deteriorates as the 
sample size increases. Moreover, the (negative) impact of increasing sample size is magnified as p 
gets smaller. 

Now consider implementing the FS algorithm for the Bernoulli mixture. Because the mixture 
has only two components, the random label switching step, y — > y', is quite simple. Indeed, we 
simply flip a fair coin. If the result is heads, then we take y' = y, and if the result is tails, then we 
take y' = y, where y denotes y with its Is and 2s flipped. The Mtm of the FS chain has entries 
given by 

k(r', s'\r, s) = - ^ vr(r', s'\y, z) iv(y\r, s, z) + - ^ vr(r', s'\y, z) n(y\r, s, z) . 



yeY ye\ 



It follows that 



k 



p m l(l-p) m / v 1 ( x X ( \ p m 0(l-p) m l / x 

p m l(1 _ p )mo Wl ^ ™2(P)+P™(l- P rm(p) w 2 (p)+p^(X- P r W0 (p) p m (l_ p )mi Wl (p) 
p^{\- p)^ Wl {p) ^(P)+P-(1-P)-^0(P) w 2 (p)+p^(X- P r W0 (p) p m„ (1 _ p) m lwi(/)) 

p m i(i-p) m o 1 s. x r \ x t \ P m o(x-p) m i 1 \ 

Note that this matrix differs from k only in the middle four elements. Indeed, the (2, 2) and (2, 3) 
elements in k have both been replaced by their average in k, and the same is true of the (3, 2) and 
(3, 3) elements. The matrix k has rank at most two, so there is at most one non-zero eigenvalue 
to find. Using the results in the Appendix along with the eigen-analysis of k performed earlier, it 
is easy to see that the non-trivial eigen-solution of k is (v±, Ai) = (i^, A2). So, the effect on the 
spectrum of adding the random label switching step is to replace the dominant eigenvalue with 0! 
(Note that Theorem Qimpl ies that A2 = Ai < X\, which justifies our ordering of the eigenvalues 
of k.) Consider again the simple numerical example with the 50:50 split of 0s and Is. In the case 
m = 10, the result of adding the extra step is to replace the dominant eigenvalue, 0.99395, by 
0.19795. When m = 20, 0.99996 is replaced by 0.15195. This suggests that, in contrast to the 
MDA algorithm, increasing sample size does not adversely affect the FS algorithm. More evidence 
for this is provided in Figure [2| which is the analogue of Figure [T] for the FS algorithm. Note that 
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the dominant eigenvalues are now substantially smaller, and no longer converge to 1 as the sample 
size increases. In fact, based on experimental evidence, it appears that, for a fixed value of p, A 2 
hits a maximum and then decreases with sample size. It is surprising that such a minor change in 
the MDA algorithm could result in such a huge improvement. In the next section, we consider a 
mixture of normal densities. 



6.2 The normal mixture 

Assume that Z\ , . . . , Z m are iid from the density 

M», r\ V ) = p V^) + (1 - P)-<P( Z -^) , 

where p G [0, 1], \i = (/xi,/i 2 ) G R 2 , r 2 = (t 2 ,t|) G R 2 , and </>(•) denotes the standard nor- 
mal density function. The prior for p is Uniform(0, 1), and the prior for (fi, r 2 ) takes the form 
7r(/ii, r 2 ) 7r(//2, t| )• As f° r 7r » we use tne standard (conditionally conjugate) prior given by 



where 7r(/ii|rf) = N(0,rf) and vr^ 2 ) = IG(2,l/2) ( |Robert and Casella| [2004| Section 9.1). 
By W ~ IG(a, 7), we mean that W is a random variable with density function proportional to 
w~ a ~ 1 exp{— 7/w}/k + (w). In contrast with the Bernoulli example from the previous subsection, 
the posterior density associated with the normal mixture is quite intractable and has a complicated 
(and uncountable) support given by X = R 2 x R 2 ^ x [0, 1]. 

The MDA algorithm is based on the complete-data posterior density, which we denote here by 
7r(/i, r 2 ,p, y\z). Again, the development in Subsection 5.1 implies that, given (fi,T 2 ,p, z), the 



elements of y are independent multinomials and the probability that the ith coordinate equals 1 
(which is one minus the probability that it equals 2) is given by 



p 



~M2 



(17) 



We sample 7r(/i, r 2 ,p\y, z) via sequential sampling from n(p\y, z) and 7r(/u, r 2 \p, y, z). The results 



in Subsection 5.1 show that p\y, z ~ Beta(ci + 1, c 2 + 1). Moreover, it's easy to show that, given 



(p, y, z), (fj>i,T 2 ) and (112, r|) are independent. Routine calculations show that 



ci + i i7 ( Cl + i; 
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and 



where z\ = ~ YaL\ I{i}(Vi) z i and s i = YT=i I{l}(Vi)( z i ~ ^i) 2 - 0f course, the distribution of 
(fj,2, rf ) given (p, y, z) has an analogous form. 

The results developed in Section [3] imply that the spectrum of the operator associated with the 
MDA chain consists of the point {0} and the eigenvalues of the Mtm of the conjugate chain, which 
lives on Y = {1, 2} m . Unfortunately, the Mtm of the conjugate chain is also intractable. Indeed, a 
generic element of this matrix has the following form: 

Hv'\y)= / / n(y'\n,T 2 ,p,z)Tr(fj 1 ,T 2 ,p\y,z)dndT 2 dp . 

JO Jr 2 + Jr 2 

This integral cannot be computed in closed form. In particular, ir(y'\fi, r 2 ,p, z) is the product of 
m probabilities of the form ( fT7] >, and the sums in the denominators of these probabilities render 
the integral intractable. However, note that k(y'\y) can be interpreted as the expected value of 
ir(y'\^L, r 2 ,p, z) with respect to the density 7r(^, r 2 ,p\y, z). Of course, for fixed z, we know how 
to draw from 7r(/i, r 2 ,p\y, z), and we have ir(y'\fi, r 2 ,p, z) in closed form. We therefore have the 
ability to estimate k(y'\y) using classical Monte Carlo. Once we have an estimate of the entire 
2 m x 2 m Mtm, we can calculate its eigenvalues. 

The same idea can be used to approximate the eigenvalues of the FS chain. The results in Sec- 
tion]?] show that we can express the FS algorithm as a DA algorithm with respect to an alternative 
complete-data posterior density, which we write as 7r*(/i, r 2 ,p, y\z). The eigenvalues of the op- 
erator defined by the FS chain are the same as those of the Mtm in which the probability of the 
transition y — > y' is given by 

/ / / TT*(y'\fi,T 2 ,p,z)Tr*(fi,T 2 ,p\y,z)dndT 2 dp . 
Jo Jr 2 + Jr 2 

It is straightforward to simulate from n*(fi, r 2 ,p\y, z), and 7r*(y'\fi, r 2 ,p, z) is available in closed 
form. 

To use our classical Monte Carlo idea to estimate the spectra associated with the MDA and FS 
chains, we must specify the data, z. Furthermore, the Bernoulli example in the previous subsection 
showed that the convergence rates of the two algorithms can depend heavily on the sample size, 
m. Thus, we would like to explore how an increasing sample size affects the convergence rates of 
the MDA and FS chains in the current context. To generate data, we simulated a random sample 
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of size 10 from a 50:50 mixture of a N(0, .55 2 ) and a N(3, .55 2 ), and this resulted in the following 
observations: 

z = Oi, . . . ,Zi ) 

= (0.2519, 2.529, -0.2930, 2.799, 3.397, 0.5596, 2.810, 2.541, 2.487, -0.1937) . 

We considered 10 different data sets ranging in size from m = 1 to m = 10. The first data set 
contained the single point z\ = 0.25192, the second contained the first two observations (zi, z%) = 
(0.25192,2.5287), the third contained (zi,z 2 ,z 3 ) = (0.25192,2.5287, -0.29303), and so on up 
to the tenth data set, which contained all ten observations. For each of these 10 data sets, we used 
the classical Monte Carlo technique described above to estimate the Mtm for both the MDA and 
FS algorithms. In particular, for each row of the Mtm we used a single Monte Carlo sample of 
size 200,000 (from ir([i, T 2 ,p|y, z) for DA, and from 7r*(^, r 2 ,p\y, z) for FS) to estimate each of 
the entries in that row. We then calculated the eigenvalues of the estimated Mtms and recorded 
the largest one. The results are shown in Figure [3] which has some interesting features. Note that 
the dominant eigenvalues of MDA chain are much closer to 1 than the corresponding dominant 
eigenvalues of the FS chain. Even at m = 5, the dominant eigenvalue of the MDA chain is already 
above 0.99. As in the previous example, the convergence rate of the MDA chain deteriorates as m 
increases. It is not clear whether the FS chain slows down as m increases. It may be the case that 
the FS eigenvalue would eventually level off, or perhaps the FS chain would eventually begin to 



speed up, as in the Bernoulli example. Note that, as proven in Subsection 5.2 when m = 1, the FS 
eigenvalue is 0. (To ascertain the accuracy of our estimates, we repeated the entire classical Monte 
Carlo simulation 6 times, with different random number seeds, and based on this, we believe that 
our eigenvalue estimates are correct up to three decimal places.) 

In the case where all 10 observations are considered, the dimension of the Mtms is 1024 x 
1024, and each element must be estimated by classical Monte Carlo. Thus, while it would be 
very interesting to consider larger sample sizes (beyond 10), and even mixtures with more than 2 
components, the matrices become quite unwieldy. 

We simulated a second set of 10 observations from the same 50:50 mixture and repeated the 
entire process for the purpose of validation. The second simulation resulted in the following data: 

z = (zi, zio) 

= (0.6699, 3.408, 0.1093, 3.289, -0.1407, 3.525, 2.454, 0.2716, -0.7443, 3.570) . 



28 



Figure [4] is the analogue of Figure [3] for the second simulation. The results are nearly identical to 
those from the first simulation. 



7 Appendix 



Consider a Mtm of the form 

abbe 
d e f f 
d f e <£ 
abbe 

and assume that all of the elements are strictly positive, so the corresponding Markov chain is 



M 



irreducible and aperiodic. Note that both of the Mtms studied in Subsection 6.1 have this form. 
Routine manipulation shows that M is reversible with respect to (tti, tt2, ^3, ni) T where tt\ = 
ad/ (ad + 2ab + cd), 112 = bni/d, tt^ = 112 and tt^ = cK\ja. In the remainder of this section, we 
perform an eigen-analysis of the matrix M. 

Of course, since M is a Mtm, it satisfies kvo = XqVq where vq = 1 and Ao = 1- Furthermore, 
since the first and fourth rows are equal, there is at least one eigenvalue equal to zero. Indeed, 
Mvs = 0, where = (c, 0, 0, — a) T . We now identify the other two eigen-solutions of M. Let 
vi = (0, 1, -1, 0) T and note that 

Mv x = (e - f)vi , 

so Ai = (e — /) is an eigenvalue. If e = /, then the middle two rows of M are equal and the rank 
of M is at most 2. (Note that Ai could be negative, implying that the operator defined by M is not 
always positive.) 

Now, let V2 = (a, 1, 1, a) , where a is a constant to be determined, and note that 

aa + 26 + ac 
ad + e + f + af 
ad + e + f + 
aa + 26 + etc 



Mv 2 



If V2 is an eigenvector with corresponding eigenvalue A2, then the first element of Mv2 must equal 
a\2\ that is, 

aa + 26 + ac = a\2 ■ 
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Now, using the fact that 26 = 1 — a — c, we have 

(a — l)(a + c) + 1 = aA 2 , 

and it follows that 

A 2 = ( "- 1)(a + c) + 1 . (18) 
a 

Again, if V2 is an eigenvector with corresponding eigenvalue A2, then the second element of Mv2 
must equal A2, or 

. , . cd 

\2 = Oid + e + j+ a — . 

a 



Now, using the fact that e = l — d — f — ^,we have 



A 2 = ^(a- l)(a + c) + 1 . 
Setting our two expressions for A 2 equal yields: 

ad(a — l)(a + c) + aa = a(a — l)(a + c) + a . 

This quadratic in a has two roots: a = 1 and 

a(a + c — 1) 



a 



d(a + c) 

The second solution is negative and corresponds to a nontrivial eigenvector. The corresponding 
eigenvalue is 

A2 = ^(o + c)(a — d) . 
If a = d, then the sum of the middle two rows of M is equal to twice the first row. 
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eigenvalues vs. sample size 
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Figure 1: The behavior of the dominant eigenvalue for the MDA chain in the Bernoulli model. The 
graph shows how the dominant eigenvalue of the MDA chain changes with sample size, m, for 
several different values of p, in the case where half the z^s are and the other half are 1. (Only even 
sample sizes are considered.) The red, blue, brown and green lines correspond to p values of 1/10, 
1/5, 1/3, and 9/20, respectively. 
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eigenvalues vs. sample size 
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Figure 2: The behavior of the dominant eigenvalue for the FS chain in the Bernoulli model. The 
graph shows how the dominant eigenvalue of the FS chain changes with sample size, m, for several 
different values of p, in the case where half the ZjS are and the other half are 1. (Only even sample 
sizes are considered.) The red, blue, brown and green lines correspond to p values of 1/10, 1/5, 
1/3, and 9/20, respectively. 



eigenvalues vs. sample size 
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Figure 3: The behavior of the dominant eigenvalue for the MDA and FS chains in the normal model. 
The graph is based on the first simulated data set and shows how the dominant eigenvalue changes 
with sample size, m, for the MDA algorithm (red line) and the FS algorithm (blue line). 
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eigenvalues vs. sample size 
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Figure 4: The behavior of the dominant eigenvalue for the MDA and FS chains in the normal model. 
The graph is based on the second simulated data set and shows how the dominant eigenvalue changes 
with sample size, m, for the MDA algorithm (red line) and the FS algorithm (blue line). 
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